function [Inew]=PriRegSENSE_LSBB(Kdata,PriI,Isens,Ksi,Par)
%This function is for the prior information regularized SENSE
%This version is only for 2D images
%INPUT
%   Kdata: partial k-space data [rnPE nFE nCH]
%   PriI: Prior information of the image, a 2D matrix, [nPE nFE ]
%   Isens: sensitivity maps [nPE nFE nCH]
%   Ksi : noise correlation matrix
%   Par: parameters
%       .priLmd: weight Lamdba for the prior informaiton, 0.2
%       .innitr: number of iternal iterations, 15
%       .eTol: determine if it is converged, 1e-4
%OUTPUT
%   Inew: reconstructed image
%-----------------------------------------------------------------
%               Feng Huang, Aug 13, 2010
%      Based on [1] Barzilai J, Borwein JM. Two-Point Step Size Gradient
%      Methods. IMA Journal of Numerical Analysis 1988;8(1):141-148
%   1. This version is trying to take advantage of initial result to
%       reduce reconstruction time
%Aug 26, 2010: The 3D version is added
%-----------------------------------------------------------------
if ndims(Kdata)==3
    %----------------------------------------------------------
    %       Preparation
    %----------------------------------------------------------
    [nrPE,nFE,nCH]=size(Kdata);
    [nPE,nFE]=size(PriI);
    %produce aliasing matrix
    [wrap_m]=sense_alias_matrix(nPE,nrPE);
    %produce the reduced FOV wrapped image
    wrapped_image = zeros(nrPE,nFE,nCH);
    for ch =1:nCH
        wrapped_image(:,:,ch) = K2Img(Kdata(:,:,ch));
        nWrImg(ch) = norm(wrapped_image(:,:,ch),'fro');
    end
    %noise correlation for whitening
    invKsi = inv(Ksi);
    %define parameters
    aP = Par.priLmd;
    innitr = Par.innitr;
    eTol = Par.eTol;
    %----------------------------------------------------------
    % Initialization
    %----------------------------------------------------------
    iter = 0;
    U = PriI;
    OBJ = realmax;
    RelErr = 1;
    CV = aP*PriI;
    bb = 1;
    N=numel(PriI);
    
    %----------------------------------------------------------
    %       Reconstruction
    %----------------------------------------------------------
    while RelErr > eTol && iter < innitr
        %calculate error, which is the fidelity term
        Err = zeros(nrPE,nFE,nCH);
        for ch=1:nCH
            %produce wrapped image using reconstruction
            tmpWrap = wrap_m*(U.*Isens(:,:,ch));
            %Energy normalization
            tmpWrap = tmpWrap/norm(tmpWrap,'fro')*nWrImg(ch);
            %calculate difference
            Err(:,:,ch) = tmpWrap - wrapped_image(:,:,ch);
        end
        %calculate norm of error
        %BB = norm(Err(:))^2/N;
        BB = norm(Err(:))^2;
        %produce FOV error image
        WrapErr = zeros(nPE,nFE,nCH);
        for PE = 1:nrPE
            wPE = find(wrap_m(PE,:));
            Rf = max(size(wPE));
            WrapErr(wPE,:,:) = repmat(Err(PE,:,:),[Rf 1 1]);
        end
        %calculate error in image space
        Grad = zeros(nPE,nFE);
        for ch=1:nCH
            Grad = Grad + WrapErr(:,:,ch).*conj(Isens(:,:,ch))*invKsi(ch,ch);
        end
        
        
        if iter > 0
            DiffU = U - Uprev;
            DiffGrad = Grad - Gradprev;
            Numer = AtimesB(DiffGrad,DiffU);
            bb = Numer/norm(DiffU,'fro')^2;
            %         fprintf('BB Step Size: %.3f\n',bb)
            if (bb < 1e-10) || (bb > 1e+10)
                fprintf('Oops, BB stepsize so small/big: %e\n',bb)
                bb = 2;
            end
        end;
        
        OBJprev = OBJ;
        %     OBJ1 = norm(U-V,'fro')^2; OBJ2=BB;
        OBJ = .5*(aP*norm(U-PriI,'fro')^2 + BB);
        RelErr = abs(OBJ-OBJprev)/OBJ;
        %fprintf('InnItr=%d, OBJ1=%.3f, OBJ2=%.3f, OBJ=%.3f\n',iter, OBJ1, OBJ2, OBJ);
        Uprev = U;
        
        % Update U
        U = (bb*U + CV - Grad)/(bb+aP);
        %     fprintf('Unorm=%.3f, Vnorm=%.3f, Grad=%.3f\n',max(abs(U(:))),max(abs(V(:))),max(abs(Grad(:))));
        
        iter = iter + 1;
        
        Gradprev = Grad;
    end
    Out.BB = BB;
    Inew = U;
elseif ndims(Kdata)==4
    %----------------------------------------------------------
    %       Preparation
    %----------------------------------------------------------
    [rnPE,rnSL,nFE,nCH]=size(Kdata);
    [nPE,nSL,nFE]=size(PriI);
    %produce aliasing matrix
    [wrap_m]=sense_alias_matrix([nPE nSL],[rnPE rnSL]);
    %produce the reduced FOV wrapped image
    wrapped_image = zeros(rnPE,rnSL,nFE,nCH);
    for ch =1:nCH
        wrapped_image(:,:,:,ch) = K2Img(Kdata(:,:,:,ch));
        tmp = wrapped_image(:,:,:,ch);
        nWrImg(ch) = norm(tmp(:),'fro');
    end
    %noise correlation for whitening
    invKsi = inv(Ksi);
    %define parameters
    aP = Par.priLmd;
    innitr = Par.innitr;
    eTol = Par.eTol;
    %----------------------------------------------------------
    % Initialization
    %----------------------------------------------------------
    iter = 0;
    U = PriI;
    OBJ = realmax;
    RelErr = 1;
    CV = aP*PriI;
    bb = 1;
    N=numel(PriI);
    
    %----------------------------------------------------------
    %       Reconstruction
    %----------------------------------------------------------
    while RelErr > eTol && iter < innitr
        %calculate error, which is the fidelity term
        Err = zeros(rnPE,rnSL,nFE,nCH);
        for ch=1:nCH
            tmpWrap = zeros(rnPE,rnSL,nFE);
            for FE = 1:nFE
                %produce wrapped image using reconstruction
                tmpWrap(:,:,FE) = wrap_m{1}*(U(:,:,FE).*Isens(:,:,FE,ch))*wrap_m{2};            
            end
            %Energy normalization
            tmpWrap = tmpWrap/norm(tmpWrap(:),'fro')*nWrImg(ch);
            %calculate difference
            Err(:,:,:,ch) = tmpWrap - wrapped_image(:,:,:,ch);
        end
        %calculate norm of error
        %BB = norm(Err(:))^2/N;
        BB = norm(Err(:))^2;
        %produce FOV error image
        WrapErr = zeros(nPE,nSL,nFE,nCH);
        for SL = 1:rnSL
            wSL = find(wrap_m{2}(:,SL));
            RfSL = max(size(wSL));
            for PE = 1:rnPE
                wPE = find(wrap_m{1}(PE,:));
                RfPE = max(size(wPE));
                WrapErr(wPE,wSL,:,:) = repmat(Err(PE,SL,:,:),[RfPE RfSL 1 1]);
            end
        end
        %calculate error in image space
        Grad = zeros(nPE,nSL,nFE);
        for ch=1:nCH
            Grad = Grad + WrapErr(:,:,:,ch).*conj(Isens(:,:,:,ch))*invKsi(ch,ch);
        end
        
        
        if iter > 0
            DiffU = U - Uprev;
            DiffGrad = Grad - Gradprev;
            Numer = AtimesB(DiffGrad,DiffU);
            bb = Numer/norm(DiffU(:),'fro')^2;
            %         fprintf('BB Step Size: %.3f\n',bb)
            if (bb < 1e-10) || (bb > 1e+10)
                fprintf('Oops, BB stepsize so small/big: %e\n',bb)
                bb = 2;
            end
        end;
        
        OBJprev = OBJ;
        %     OBJ1 = norm(U-V,'fro')^2; OBJ2=BB;
        OBJ = .5*(aP*norm(U(:)-PriI(:),'fro')^2 + BB);
        RelErr = abs(OBJ-OBJprev)/OBJ;
        %fprintf('InnItr=%d, OBJ1=%.3f, OBJ2=%.3f, OBJ=%.3f\n',iter, OBJ1, OBJ2, OBJ);
        Uprev = U;
        
        % Update U
        U = (bb*U + CV - Grad)/(bb+aP);
        %     fprintf('Unorm=%.3f, Vnorm=%.3f, Grad=%.3f\n',max(abs(U(:))),max(abs(V(:))),max(abs(Grad(:))));
        
        iter = iter + 1;
        
        Gradprev = Grad;
    end
    Out.BB = BB;
    Inew = U;
end
return;